function stats = calibrate_skill_dist_post(version, varargin)
% post-calibration analysis - generates files with conditional
% densities and changes in wage-distribution deciles. Plotting is
% done only for immediate diagnostics. Plots for paper are produced by
% R scripts.

addpath('../Densities')
addpath('../UtilityFunctions')
addpath('../../matlab')

dir_in = '../../../data/generated/matlab/Calibration/';

load([dir_in,version,'_results.mat'])
par_hat = sol;

valid_num_nodes = @(x) isnumeric(x) && x==floor(x) && (x>0);

p = inputParser;
p.addOptional('lb',[],@isstruct)
p.addOptional('ub',[],@isstruct)
p.addParameter('plot_partials',false,@islogical)
p.addParameter('num_nodes_f',cal.num_nodes_f,valid_num_nodes)
p.addParameter('num_nodes_g',cal.num_nodes_g,valid_num_nodes)
p.addParameter('num_nodes_bins',cal.num_nodes_bins,valid_num_nodes)
p.addParameter('num_nodes_phi',cal.num_nodes_phi,valid_num_nodes)

p.parse(varargin{:})
lb = p.Results.lb;
ub = p.Results.ub;
num_nodes_f = p.Results.num_nodes_f;
num_nodes_g = p.Results.num_nodes_g;
num_nodes_bins = p.Results.num_nodes_bins;
num_nodes_phi = p.Results.num_nodes_phi;
plot_partials = p.Results.plot_partials;

w_lb = cal.w_lb;
w_ub = cal.w_ub;

dist_0 = cal.dist;
dist_1 = cal.dist;

stats = cal.moments_skill_prices(par_hat);

figure(1)
subplot(1,3,1)
bar([cal.targets.freqM.val,...
    dist_0.h_M_mass_wagebins(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0, ...
    cal.wagebins.M, stats.participation_rate_0,par_hat)])
ylim([0,0.1])
title('M')

subplot(1,3,2)
bar([cal.targets.freqR.val,...
    dist_0.h_R_mass_wagebins(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0, ...
    cal.wagebins.R, stats.participation_rate_0,par_hat)])
ylim([0,0.1])
title('R')

subplot(1,3,3)
bar([cal.targets.freqC.val,...
    dist_0.h_C_mass_wagebins(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0, ...
    cal.wagebins.C, stats.participation_rate_0,par_hat)])
ylim([0,0.1])
title('C')
legend('data','model')


%% summary statistics
stats.emp_M_0 = dist_0.emp_M(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,w_lb,w_ub,stats.participation_rate_0,par_hat);
stats.emp_R_0 = dist_0.emp_R(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,w_lb,w_ub,stats.participation_rate_0,par_hat);
stats.emp_C_0 = dist_0.emp_C(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,w_lb,w_ub,stats.participation_rate_0,par_hat);

stats.emp_M_1 = dist_1.emp_M(par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,w_lb,w_ub,stats.participation_rate_1,par_hat);
stats.emp_R_1 = dist_1.emp_R(par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,w_lb,w_ub,stats.participation_rate_1,par_hat);
stats.emp_C_1 = dist_1.emp_C(par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,w_lb,w_ub,stats.participation_rate_1,par_hat);

handle_f_0 = @(w) dist_0.dens_skill.f(w,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat);
stats.total_mass = tools.Integration.integral_legendre(handle_f_0,w_lb, ...
    w_ub,cal.num_nodes_f);

stats.avg_wage_M_0 = dist_0.average_wage_M(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,w_lb,w_ub,stats.participation_rate_0,par_hat);
stats.avg_wage_R_0 = dist_0.average_wage_R(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,w_lb,w_ub,stats.participation_rate_0,par_hat);
stats.avg_wage_C_0 = dist_0.average_wage_C(par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,w_lb,w_ub,stats.participation_rate_0,par_hat);


stats.avg_wage_M_1 = dist_1.average_wage_M(par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,w_lb,w_ub,stats.participation_rate_1,par_hat);
stats.avg_wage_R_1 = dist_1.average_wage_R(par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,w_lb,w_ub,stats.participation_rate_1,par_hat);
stats.avg_wage_C_1 = dist_1.average_wage_C(par_hat.Y_M_1, ...
    par_hat.Y_R_1,par_hat.Y_C_1,w_lb,w_ub,stats.participation_rate_1,par_hat);

% percentiles
handle_h_0 = @(w) 1/stats.participation_rate_0 .* dist_0.h(w,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat);

pct = 0.1:0.1:0.9;
fprintf('\nusing interpolation\n')
tic
[stats.percentiles_0, nodes, weights] = tools.Integration.compute_percentiles(handle_h_0,pct,1,200,100);
toc


handle_h_1 = @(w) 1/stats.participation_rate_1 .* dist_1.h(w,par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,par_hat);

fprintf('\nusing interpolation\n')
tic
stats.percentiles_1 = tools.Integration.compute_percentiles(handle_h_1,pct,1,200,100,nodes,weights);
toc

stats.percentiles_change = stats.percentiles_1 - stats.percentiles_0;
figure(3)
bar(stats.percentiles_change)
title('Change in model deciles between periods')

% percentile change by sector
% idea is to construct distribution at t=1 step by step, starting
% R, then M, then C. First construct handles at 0 and 1 by occ

handle_h_0_M = @(w) 1/stats.participation_rate_0 .* dist_0.h_M(w,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat);
handle_h_0_R = @(w) 1/stats.participation_rate_0 .* dist_0.h_R(w,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat);
handle_h_0_C = @(w) 1/stats.participation_rate_0 .* dist_0.h_C(w,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat);

handle_h_1_M = @(w) 1/stats.participation_rate_1 .* dist_1.h_M(w,par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,par_hat);
handle_h_1_R = @(w) 1/stats.participation_rate_1 .* dist_1.h_R(w,par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,par_hat);
handle_h_1_C = @(w) 1/stats.participation_rate_1 .* dist_1.h_C(w,par_hat.Y_M_1,par_hat.Y_R_1,par_hat.Y_C_1,par_hat);

% now, mix in changes by occupation step by step
handle_h_1_R_only = @(w) handle_h_1_R(w) + handle_h_0_M(w) + handle_h_0_C(w);
handle_h_1_R_M = @(w) handle_h_1_R(w) + handle_h_1_M(w) + handle_h_0_C(w);
handle_h_1_R_M_C = @(w) handle_h_1_R(w) + handle_h_1_M(w) + handle_h_1_C(w);

% compute percentiles
stats.percentiles_1_R_only = tools.Integration.compute_percentiles(handle_h_1_R_only,pct,1,200,100,nodes,weights);
stats.percentiles_1_R_M = tools.Integration.compute_percentiles(handle_h_1_R_M,pct,1,200,100,nodes,weights);
stats.percentiles_1_R_M_C = tools.Integration.compute_percentiles(handle_h_1_R_M_C,pct,1,200,100,nodes,weights);

% compute percentile changes
stats.percentiles_change_R_only = stats.percentiles_1_R_only - stats.percentiles_0;
stats.percentiles_change_R_M = stats.percentiles_1_R_M - stats.percentiles_0;
stats.percentiles_change_R_M_C = stats.percentiles_1_R_M_C - stats.percentiles_0;

%% plot against wage data histograms

cps_data = tools.csv2struct(['../../..//data/' ...
    'generated/r/morg_data_stacked_90.csv']);


tic
w_range = 1:0.01:200;
figure(5)
subplot(1,3,1)
histogram(cps_data.ln_w_pareto(cps_data.ind1==1&cps_data.wagesample==1),'Normalization','pdf')

dat.ln_w = log(w_range)';
dat.cond_dens_M = (w_range.*1/stats.participation_rate_0 .*dist_0.h_M(w_range,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat)./stats.emp_M_0)';
dat.cond_dens_R = (w_range.*1/stats.participation_rate_0 .*dist_0.h_R(w_range,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat)./stats.emp_R_0)';
dat.cond_dens_C = (w_range.*1/stats.participation_rate_0 .*dist_0.h_C(w_range,par_hat.Y_M_0,par_hat.Y_R_0,par_hat.Y_C_0,par_hat)./stats.emp_C_0)';

title('occ 1')
hold on
plot(log(w_range),dat.cond_dens_M,'LineWidth',3)
hold off

subplot(1,3,2)
histogram(cps_data.ln_w_pareto(cps_data.ind2==1&cps_data.wagesample==1),'Normalization','pdf')
title('occ 2')
hold on
plot(log(w_range),dat.cond_dens_R,'LineWidth',3)
hold off

subplot(1,3,3)
histogram(cps_data.ln_w_pareto(cps_data.ind3==1&cps_data.wagesample==1),'Normalization','pdf')
title('occ 3')
hold on
plot(log(w_range),dat.cond_dens_C,'LineWidth',3)
hold off

tools.struct2csv(dat,['../../../data/generated/matlab/Calibration/',version,'_conditional_densities.csv']);

toc


%% partial derivatives
if plot_partials
    opt = cal.objective_skill_dist;
    Diag = DiagnosticsForOptimizationProblems(opt, lb, ub);
    Diag.plot_objective_all_variables(par_hat, 20)
end

%% compare wage-decile elasticities in model and data

figure(6)
bar([cal.targets.eps_w_K.val,stats.percentiles_change./ ...
    stats.percentiles_0./stats.K_B_rel_change])
legend('data','model')

deciles_eps.data = cal.targets.eps_w_K.val;
deciles_eps.model = stats.percentiles_change./ ...
    stats.percentiles_0./stats.K_B_rel_change;

tools.struct2csv(deciles_eps,['../../../data/generated/matlab/Calibration/',version,'_deciles_eps.csv']);

%% plot heathcote tax function
mtax_y = @(y) cal.dist.mtax_hsv_y(y, par_hat.lambda_hsv, ...
    par_hat.tau_hsv);
eps = cal.dist.util.epsilon;
y_w = @(w) cal.dist.scale_inc.*cal.dist.y_hsv(w, eps, par_hat.lambda_hsv, ...
    par_hat.tau_hsv, par_hat.xi);

wrange = 1:100;
y = y_w(wrange);
y_scaled = y./cal.dist.scale_inc./stats.avg_labor_income;
mtax = mtax_y(y);

figure(7)
plot(y_scaled,mtax)
title('Marginal Tax Rates, x=1 is avg income')

%% decompose percentile changes
ylim_min = min(min(stats.percentiles_change_R_only./ ...
    stats.percentiles_0./stats.K_B_rel_change, ...
               cal.targets.eps_w_K.val));

ylim_max = max(max(stats.percentiles_change_R_only./ ...
    stats.percentiles_0./stats.K_B_rel_change, ...
               cal.targets.eps_w_K.val));

figure(7)
bar([cal.targets.eps_w_K.val,stats.percentiles_change_R_only./ ...
    stats.percentiles_0./stats.K_B_rel_change])
legend('data','model')
title('change due to R')

figure(8)
bar([cal.targets.eps_w_K.val,stats.percentiles_change_R_M./ ...
    stats.percentiles_0./stats.K_B_rel_change])
legend('data','model')
title('change due to R and M')
ylim([ylim_min, ylim_max])

figure(9)
bar([cal.targets.eps_w_K.val,stats.percentiles_change_R_M_C./ ...
    stats.percentiles_0./stats.K_B_rel_change])
legend('data','model')
title('change due to R, M, C')
ylim([ylim_min, ylim_max])

figure(10)
bar(stats.percentiles_change_R_only)
title('Change in model deciles between periods, R only')


end
